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The emergence of the bipolaronic phase and the formation of the heavy-electron state in the anharmonic Holstein 
model are investigated using the dynamical mean-field theory in combination with the exact diagonalization method. 
For a weak anharmonicity, it is confirmed that the first-order polaron-bipolaron transition occurs from the observation 
of a discontinuity in the behavior of several physical quantities. When the anharmonicity is gradually increased, the 
polaron-bipolaron transition temperature is reduced as well as the critical values of the electron-phonon coupling con- 
stant for polaron-bipolaron transition. For a strong anharmonicity, the polaron-bipolaron transition eventually changes to 
a crossover behavior. The effect of anharmonicity on the formation of the heavy-electron state near the polaron-bipolaron 
transition and the crossover region is discussed in detail. 
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1. Introduction 

Recently, heavy-electron phenomena have attracted re- 
newed attention in the research field of condensed matter 
physics. 1 A traditional mechanism of the emergence of the 
heavy-electron state is based on quantum criticality induced 
by the competition between the Kondo effect and Ruderman- 
Kittel-Kasuya-Yosida interaction. The Kondo effect due to 
local magnetic moment has been well understood, 2 but the 
Kondo-like phenomenon occurs in a more general manner, 
when a localized entity with internal degrees of freedom is 
embedded in a conduction electron system and quantum- 
mechanical exchange interaction effectively works between 
local degrees of freedom and conduction electrons. 

In particular, the Kondo phenomenon with a phonon ori- 
gin has been potentially discussed. First, Kondo himself has 
considered a conduction electron system coupled with a lo- 
cal double-well potential. 3,4 Having two possible electron po- 
sition in a double-well potential play roles in pseudo-spins, 
and the Kondo-like behavior is considered to appear in such 
a two-level system. In fact, it has been shown that the two- 
level Kondo system exhibits the same behavior as the mag- 
netic Kondo effect. 5 7 Another important issue regarding the 
Kondo effect with a phonon origin has been shown by Yu and 
Anderson. 8 They have pointed out that the scattering process 
between spinless s-wave conduction electrons to p-wave ones 
is induced by ion displacement. The model proposed by Yu 
and Anderson has been shown to be mapped to the two-level 
Kondo model at low temperatures. 9 ' 10 

A recent revival of research on the Kondo effect with a 
phonon origin has been triggered by active experimental in- 
vestigations on cage structure compounds such as filled skut- 
terudites, 11-14 clathrates, 15-20 and /3-pyrochlore oxides. 21-25 
In these materials, a guest atom is contained in a cage com- 
posed of relatively light atoms and oscillates with a large am- 
plitude in a potential with a strong anharmonicity. Such local 
oscillation with a large amplitude is called rattling, and exotic 
phenomena induced by rattling have attracted much attention 
in the research of strongly correlated electron materials with 



a cage structure. 

An example of such active investigations is the study of 
a magnetically robust heavy-electron behavior observed in 
SmOs4Sbi2. 14 The electronic specific heat coefficient j c is 
in proportion to the effective mass of electrons, but in this 
material, 7 e is several hundred times larger than that of a free- 
electron system and is found to be almost unchanged even if 
a magnetic field up to 30 Tesla is applied. If the heavy ef- 
fective mass originates from the quantum criticality due to 
the traditional Kondo effect with a magnetic origin, -f c should 
be significantly suppressed by the applied magnetic field. The 
origin of the heavy-electron state in SmOs4Sbi2 is considered 
to be electron-rattling interaction. In fact, four- and six-level 
Kondo systems have been analyzed and the magnetically ro- 
bust heavy-electron state has been actually obtained. 26-28 The 
periodic Anderson-Holstein model has been analyzed using 
the dynamical mean-field theory, and the mass enhancement 
due to large lattice fluctuations and phonon softening towards 
a double-well potential have been addressed. 29-31 

Kondo phenomena in the conduction electron system 
coupled with local Jahn-Teller phonons 32 ' 33 and Holstein 
phonons 34 ' 35 have been discussed for the realization of a non- 
magnetic Kondo effect. From the numerical evaluation of 7 e 
in the conduction electron system coupled with local anhar- 
monic Holstein phonons, it has been shown that 7 C is largely 
enhanced by rattling and is actually robust under an applied 
magnetic field. 36 Furthermore, it has been pointed out that the 
Kondo effect due to rattling should exhibit a peculiar isotope 
effect, which is experimental evidence of rattling-induced 
heavy-electron phenomena. 37 Quite recently, the vibration of 
magnetic ions has been explicitly included in the spinful Yu- 
Anderson model, and the two-channel Kondo effect has been 
comprehensively discussed. 38 ^ 10 

When we turn our attention to the /3-pyrochlore oxide 
KOS2O6, significant effects of the rattling of a K ion in a tetra- 
hedral cage composed of O and Os ions have been discussed. 
The rattling-associated anomaly is found in the form of the 



first-order structural transition at T, 



7.5 K, which is dif- 
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ficult to be influenced by a magnetic field and is not accom- 
panied by any symmetry breaking. 23 ' 24 It has been reported 
that electric resistivity is in proportion to T 2 for T < T p , 
while it behaves as oc >/T for T > T p . In the peculiar be- 
havior at high temperatures, it has been suggested that anhar- 
monic phonons play important roles. 41 The transition at T p 
has been discussed in the context of electron-rattling inter- 
action. A possible liquid-gas-type rattling transition and the 
multipole transition driven by the octupole component of K 
ion rattling have been pointed out. 42 ' 43 From the analysis of 
the harmonic Holstein model, it has been reported that the 
first-order transition from the polaronic state to the bipola- 
ronic state occurs in the strong-interaction region at low tem- 
peratures. 44 ^ 7 

In this study, to obtain deep insight into the effect of anhar- 
monicity on the polaron-bipolaron transition and the heavy- 
electron state, we analyze the anharmonic Holstein model us- 
ing the dynamical mean-field theory (DMFT) in combination 
with the exact diagonalization method at low temperatures 
such as 10 -4 of the conduction electron bandwidth. 48 When 
the anharmonicity is weak, i.e., for nearly harmonic phonons, 
we again find the first-order polaron-bipolaron transition due 
to the observation of a discontinuity in the behavior of phys- 
ical quantities. When the anharmonicity is increased, the 
polaron-bipolaron transition temperature is reduced and the 
critical value of electron-phonon interaction becomes smaller. 
For a strong anharmonicity, the polaron-bipolaron transition 
disappears and it becomes to a crossover behavior. We dis- 
cuss the effect of anharmonicity on the heavy-electron state 
near the polaron-bipolaron transition and crossover regions. 

The organization of this paper is as follows. In § 2, we in- 
troduce the anharmonic Holstein model and explain the an- 
harmonic potential used in this paper. We also provide a brief 
explanation of the DMFT. In § 3, we discuss the properties of 
the system obtained using the DMFT. Several physical quan- 
tities are discussed when we change the electron-phonon in- 
teraction, the anharmonicity of phonons, and temperature. Fi- 
nally, in § 4, we summarize this paper. Throughout this paper, 
we use such units as H = ks = 1. 



2. Model and Method 

In this section, we introduce the Hamiltonian as 



H — £ fe c fecr C feCT + -ffeph, 



(1) 



where is the energy of conduction electrons, c£ ff is a cre- 
ation operator for conduction electrons with a wave vector 
k and a spin a, and -ff ep h denotes a local electron-phonon 
term. In the following, we describe H ep h and the potential for 
the vibration of the guest atom. Then, we briefly explain the 
scheme of the DMFT to solve the Hamiltonian. 

2.1 Electron-phonon coupling term 

We consider a situation in which electrons are coupled with 
the local oscillation of an atom. Such a situation is expressed 
by 



H, 



eph 



E 



9Qi{rH-l) + ^+V(Qi) 



(2) 



where g denotes the coupling constant between electron den- 
sity and the oscillation of the atom, i indicates the atomic site, 
n i = 2~2 a c L Cicr ' CiCT 1S tne annihilation operator of an elec- 
tron with a spin a at a site i, Qi is the normal coordinate of 
breathing-mode oscillation of the atom, Pi indicates the corre- 
sponding canonical momentum, M is mass of oscillator atom, 
and V{Qi) is a potential for atom, which is expressed by 



V(Qi) = k 2 Q 2 + k 4 Qt- 



(3) 



Here, k 2 and fc 4 respectively denote the coefficients for the 
second- and fourth-order terms of the potential for the atom. 

Note that in the coupling between electrons and oscillation, 
we subtract unity from the electron number for convenience 
in the numerical calculation. If we do not carry out this sub- 
traction, the electron-phonon coupling is effectively enhanced 
at doubly occupied sites and it is necessary to prepare a large 
number of phonon basis to obtain convergent results. Since we 
are interested in the bipolaronic state, it is crucial to perform 
the calculation with a significant amount of double occupancy 
with sufficient precision. Thus, we use the electron-phonon 
term in eq. (2). 

For more calculations, it is convenient to introduce the 
second-quantized phonon operator by following the stan- 
dard quantum mechanical procedure. The displacement is ex- 
pressed by Qi = Qoipi + b\), where Qo = 1/\ / 2Mujq, ujq 
is the phonon energy, and hi is the annihilation operator of 
phonons. Note that Qo denotes the amplitude of zero-point 
oscillation for harmonic phonons. Since we consider the an- 
harmonic oscillation including the case of a negative k 2 , we 
do not impose an explicit relation between k 2 and uj - 

Then, we rewrite i7 cp h using phonon operators as 



H C ph — ^0 



[yS(&i + 6t)(n i -l) + &t& i + l/2 



(4) 



+ /32(b i + b\) 2 + (3 i (b i + bl) 4 ], 



where a, 2 , and 04 are given by 

2k 2 



a ■■ 
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02 



(5) 



Among them, a indicates the nondimensional electron- 
phonon coupling constant, while 2 and 04 denote the nondi- 
mensional second- and fourth-order anharmonicity parame- 
ters, respectively. Note that the anharmonicity of the potential 
is controlled by adjusting 2 , which becomes both positive 
and negative, while we always set 04 > to restrict the os- 
cillation of an atom in a finite space. We also note that the 
harmonic case is denoted by 2 = 4 = 0. 

2.2 Anharmonic potential 

Now, we discuss the anharmonic potential for the oscilla- 
tion of an atom. For this purpose, it is convenient to intro- 
duce the nondimensional displacement qi through the relation 
Qi = Qi/ Qo in accordance with the second-quantization of 
a phonon. Using nondimensional displacement, we obtain the 
potential as 

wo ( 02 + ] ) q- +04qf 



V{qi 



(6) 



Note that the energy scale of V(qi) is taken as ujq. 
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Fig. 1. (Color online) (a) Local phonon potential V(q) for several values 
of /32 from to —0.7 with 04 = 0.02. (b)-(e) Several typical potentials with 
02 = 0, —0.25, —0.5, and —0.7, respectively. We show eigenenergy levels 
with dotted horizontal lines. 



(P 4 =0.02) 



first excitation energy 
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third excitation energy 
potential well depth 




Fig. 2. (Color online) Energies in the local phonon potentials. The first, 
second, and third excitation energies of the local phonon and the depth of the 
potential well are shown. The phonon frequency ujq = 0.4 is also indicated. 



crease /3 2 , the first excitation energy is gradually decreased 
and eventually at /3 2 ~ —0.7, it becomes almost zero. On the 
other hand, the second and third excitation energies are rather 
increased in the region of fa < —0.5, since the ground-state 
energy is rapidly decreased. 



In Fig. 1(a), we show several anharmonic potentials by 
changing /3 2 for /3 4 = 0.02. For /3 2 > -0.25, the poten- 
tial minimum is always situated at qi =0 and the shape of 
a single-well potential is similar to that for a harmonic po- 
tential. However, in the region of > /3 2 > —0.25, the 
bottom of the potential is found to be wide and flat, since 
at /3 2 = —0.25, the second-order term of the potential van- 
ishes. When /3 2 is decreased from —0.25, the potential shape 
is changed, since potential minima appear at ±g m ; n , given by 
q min = \/Wz + l/4)/(2/3 4 ). For < -0.25, the potential 
minimum rapidly decrease, leading to the double-well poten- 
tial. We consider that the decrease in /3 2 indicates the increase 
in anharmonicity. 

In Figs. 1(b)- 1(e), we show some typical shapes for a 
single-well type for /3 2 = 0, a flat-bottom type for /3 2 = 
—0.25, a shallow double-well type for /3 2 = —0.5, and a deep 
double- well type for /? 2 = —0.7. We also show the eigenen- 
ergy levels obtained by the diagonalization of i/ cp h- It is ob- 
served that the energy levels tend to move to the lower side 
and that the width between adjacent levels becomes smaller, 
when /3 2 is decreased. In the double-well type cases shown in 
Figs. 1(d) and 1(e), several levels are found inside the wells. 
When we further decrease /3 2 , the levels of the ground and 
first-excited states exist deep inside the wells and the width 
between adjacent levels is extremely small. We expect that 
the ground and first-excited states finally become degenerate 
for a large negative /3 2 with a deep double-well type potential. 

In Fig. 2, we show several excitation energies and the depth 
of the potential well, defined by ^(0) — V(q m i n ). We also de- 
pict ljq = 0.4 line in the graph, which corresponds to the 
first excitation energy in the harmonic phonon case. When 
we decrease /3 2 , i.e., increase anharmonicity, we find that the 
excitation energies are totally suppressed. For /3 2 < —0.25, 
the potential wells are formed at qi — ±q m i n . When we de- 



2.3 Dynamical mean-field theory 

To solve the model Hamiltonian, it is necessary to use 
appropriate approximation depending on the problem. Here, 
we adopt the dynamical mean-field theory (DMFT). 49 In the 
present case, the model is mapped onto an effective impurity 
Anderson-Holstein model. 50 In the following, we briefly ex- 
plain the scheme of the DMFT. 

The local electron Green's function G(iui n ) is given by 

G(iu n ) = J2- : — ~ (7) 

k 



i k - A* - S(zw„) ' 



where u n is the fermion Matsubara frequency, given by cj n = 
7rT(2n + 1) with an integer n, T is the temperature, /i is the 
chemical potential, and E(iw n ) is the electron self-energy. 
Note that, in the DMFT, the momentum dependence of the 
self-energy can be ignored. 

To perform momentum summation, it is necessary to spec- 
ify the lattice type. Here, we consider the Bethe lattice with a 
semielliptic density of states (DOS) given by 

P( £ ) = -^\/W 2 -4 e 2 , (8) 

where W is the bandwidth. Using this DOS, we obtain the 
condition for the local Green's function as 



G {iu n ) 1 = iw n - n - {W/A) 2 G{ioj n ), 



(9) 



where Go(iui n ) is the bare local electron Green's function. 
In the present work, we determine Go by the effective im- 
purity Anderson-Holstein model with a — in an effective 
medium determined self-consistently. The effective impurity 
Anderson-Holstein model with a finite a is solved by the ex- 
act diagonalization method for a finite-sized cluster to obtain 



3 



J. Phys. Soc. Jpn. 



DRAFT 




ues with /?4 = 0.02 and T = 0.0025. Fig. 4. (Color online) Double occupancy d vs a for several /?2 values with 

04, = 0.02 and T = 0.0025. 



G(iuj n ) at a finite temperature, given by 

l — e -Et/T + e -E j /T 

G(iu n ) = 7 E ■ , t F- 0' , (10) 

3,1 

where Z is the partition function given by Z — Y^j e~ Ej / T ', 
\j) is the eigenstate of the effective impurity Anderson- 
Holstein model, and Ej is the corresponding eigenenergy. 

In the present paper, we use the five-site cluster and the cut- 
off number N c for the phonon basis is set to be 12. We have 
checked the convergence of the calculations in comparison 
with the results for the six-site cluster and N c — 15. We con- 
centrate our attention to the half-filling case with (m) = 1. 
We set W = 4 and cj = 0.4 in the following calculations. 
Note that we restrict ourselves to the case with the normal 
state without any symmetry breaking. 

3. Results of DMFT Calculations 

3.1 Lattice fluctuations 

Let us start our discussion about the anharmonicity depen- 
dences of various physical quantities for a fixed temperature. 
First, we discuss the anharmonicity effect on the lattice fluctu- 
ation yj {q 2 ) . In Fig. 3, we show the a dependence of the lat- 
tice fluctuation yj (q 2 ) for several values of fa with fa = 0.02 
and T = 0.0025. It is found that the curves for \J (q 2 ) are 
monotonically increasing functions. In the following, we dis- 
cuss the change in the curves for different values of fa. 

When we decrease fa, i.e., strengthen the anharmonicity, 
the amplitude of the lattice fluctuations becomes larger. In 
fact, at a = 5, the lattice fluctuation increases from 2.2 to 
~ 4.1 when fa decreases from to —0.7. Note that the lattice 
fluctuation at a = already exhibits a similar tendency as it 
increases from ~ 0.9 to ~ 3.2 when fa is changed from 
to —0.7. This tendency can be understood from the potential 
shape, as shown in Fig. 1. Namely, when fa is decreased, the 
total width of the potential becomes larger and the amplitude 
of the oscillation in the potential increases owing to thermally 
activated and quantum tunneling motions. 

Next, we comment on the change in the function shape 
due to fa. For small absolute values of fa, the shape of 
the curve is convex-downward for a small a, while it be- 



comes slightly convex-upward for a large a. For instance, for 
fa = 0, the change in the behavior can be observed at a ~ 4. 
When the absolute value of fa becomes larger, the curve be- 
comes steeper and a at which the convex-downward behavior 
changes to a slightly convex-upward one becomes smaller. Fi- 
nally, for fa — —0.6 or —0.7, the behavior for a small a also 
becomes a linear function. These changes are thought to be 
caused by the enhancement of the effective electron-phonon 
interaction coupled with the anharmonic potential. 

For a negative fa with small absolute value, we find the 
hysteresis region at approximately a where the increasing be- 
havior is changed. In actual calculations, we obtain two solu- 
tions for the same a by gradually increasing and decreasing 
a. One of the solutions is stable and the other is metastable. 
Thus, this region can be accounted for by the hysteresis be- 
havior observed in the experimental measurements as well as 
by the first-order phase transition point where the stable and 
metastable solutions are switched. 

With decreasing fa, the hysteresis region becomes smaller. 
For fa < —0.5, it eventually disappears and becomes a 
smooth crossover between the two solutions. Since the size 
of the hysteresis region is considered to be related to the en- 
ergy scale of the first-order phase transition temperature, this 
result indicates that the anharmonicity suppresses the first- 
order phase transition. Note that in the present calculations, 
the systematic reduction in the size of the hysteresis region is 
obtained at a certain accuracy. When we attempt to quantita- 
tively obtain improved results, it is necessary to resort to an 
extrapolation technique to estimate the size of the hysteresis 
region. 

3.2 Double occupancy 

Now we move to the discussion of the electronic states. 
For this purpose, we examine the behavior of the double oc- 
cupancy d = (riffii) as one of the typical physical quanti- 
ties. In Fig. 4, d is plotted as a function of a for several fa 
values with fa = 0.02 and T = 0.0025. For a = 0, we 
find d — 0.25, since four local electron states are degen- 
erate in the noninteracting case. When we increase a, d in- 
creases monotonically, but the increasing behavior changes 
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Xc 400 



200 




Fig. 5. (Color online) Local spin susceptibility Xs loc vs a f° r several ft^ 
values with /3 4 = 0.02 and T = 0.0025. 



Fig. 6. (Color online) Local charge susceptibility Xc loc vs a for several 
@2 values with /3 4 = 0.02 and T = 0.0025. 



during the transition or the crossover. For a small a, d is a 
convex-downward (almost linear) function at a large (small) 
/?2, while it is a convex-upward function for a large a. For 
a sufficiently large a, d asymptotically approaches 0.5, since 
the vacant and the double-occupied states are degenerate at 
half-filling in the large a limit. 

We remark that d directly indicates the amount of the bipo- 
laronic state in which electrons with antiparallel spins are cou- 
pled due to the effective attractive interaction mediated by 
phonons. Thus, we characterize the two solutions in the hys- 
teresis region using d. Namely, among the two solutions, we 
define the solution with a larger d as the bipolaronic solution, 
while that with a smaller d is called the polaronic solution. 
Note that, in the definition, the polaronic solution does not 
indicate d = 0. 

When we consider the f3 2 dependence, we again find sev- 
eral remarks that we have already described in the previous 
subsection. Namely, the coexistence region disappears for a 
small (3 2 and the polaron-bipolaron transition or crossover 
point decreases with decreasing f3 2 . For a small f3 2 , a steep 
increase in d with increasing a is also observed. 

3.3 Local spin and charge susceptibilities 

In the bipolaronic state with a large d, it is expected that the 
local charge fluctuation is enhanced together with the local 
lattice fluctuation, while the local moment is suppressed. In 
fact, from the definition of d in the previous subsection, we 
obtain ((n, - {n.1)) 2 ) = 2d and (s zi 2 ) = (1 - 2ef)/4, where 
s zi — ( c ]-f c it — c \i c ii)l^- Th en > we examine the local charge 
and spin susceptibilities Xc° c and Xc° c > defined as 



X 



loc 



dT(ni(T)ni), 



and 



Joe 



1/T 



dT(s zi (T)s zi }, 



(ID 



(12) 



respectively. 

In Figs. 5 and 6, we show Xs loc and Xc loc > respectively, 
for several anharmonic potentials with T = 0.0025. With 



an increasing function. The degree of decrease or increase 
tendency is enhanced, when f3 2 is decreased, i.e., the anhar- 
monicity is increased. The behavior agrees well with our ex- 
pectation of the local charge and spin fluctuations, mentioned 
above. For the solutions with a sufficiently large a, irrespec- 
tive of /3 2 , Xs loc asymptotically approaches zero, while Xc loc 
approaches a finite specific value, ~ 400. The behavior has 
been obtained in the calculations for a harmonic phonon. 47 
Thus, the present result confirms that the states with a large 
a describe the bipolaronic state, in which charge fluctuations 
are dominant, while spin fluctuations are totally suppressed. 

3.4 Reno rmalization factor 

In this subsection, we discuss the possibility of the heavy- 
electron state due to the strong electron-phonon interaction 
with the anharmonicity. For this purpose, we estimate the 
renormalization factor z defined by 

dRe£(o;) 



= 1 - 



dui 



(13) 



increasing a, Xs loc is a decreasing function, while x c loc is 



where is the electron self-energy with a real frequency 
lj. In general, z means the renormalization effect of the con- 
duction electron state. 

In Fig. 7, we plot z as a function of a for several j3 2 values 
at a low temperature T = 0.0025. With increasing a, irrespec- 
tive of /?2 values, z monotonically decreases and finally goes 
to zero. As well as other physical quantities discussed in the 
previous subsections, for a weak-anharmonicity cases with 
/?2 > —0.3, we find the hysteresis region, in which large- and 
small-z solutions coexist at the same a. On the other hand, for 
strong-anharmonicity cases with f3 2 < —0.4, the crossover 
between the two solutions is observed. The behavior of the 
decreasing z seems to depend on (3 2 - For the harmonic-like 
potential with f3 2 > —0.25, z is a slightly convex-upward 
function for a small a, while for the double-well potential 
with /?2 < —0.25, it changes to a convex-downward func- 
tion. For a smaller f3 2 , we find a steeper decrease in z for a 
small a with a polaronic solution. 

Here, we recall that the inverse of z indicates the mass en- 
hancement of electrons, expressed by = m*/m, where 
m* denotes the effective mass of the conduction electron and 
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Fig. 7. (Color online) Renormalization factor zvsa for several /?2 values 
with /3 4 = 0.02 and T = 0.0025. 



m is the bare electron mass. Thus, the small-z solution in- 
dicates the formation of a heavy-electron state. Note that the 
solution in the limit of z — >• indicates the bipolaronic state 
with a large a. The heavy-electron state should appear in the 
polaronic phase in the vicinity of the bipolaronic phase. In 
this sense, the case in which the solutions with a small z are 
widely distributed in the parameter region seems to be advan- 
tageous for the emergence of the heavy-electron state. 

For weak-anharmonicity cases such as fa — 0, small- 
z solutions with a possible heavy-electron state are found 
only in the vicinity of the two-solution coexistence region. 
Thus, there are small possibilities of the heavy-electron 
state for weak-anharmonicity cases. On the other hand, for 
strong-anharmonicity case with crossover solutions such as 
fa = —0.5, z exhibits a convex-downward behavior in the 
crossover region. Namely, small-z solutions are distributed 
in the relatively wide region of a. When fa is decreased, 
the tail-like structure can be found in the z behavior in the 
crossover region, indicating that the region of small-z solu- 
tions becomes wide. Therefore, the heavy-electron state due 
to the electron-phonon interaction is brought about by the 
large anharmonicity. 

To discuss the heavy-electron state in more detail, we show 
the phase diagram concerning z in the (a, fa) phase for 
T = 0.0025 in Fig. 8. Here, we define that the heavy-electron 
state is characterized by the solution with 10 < z _1 < 1000 
and that the corresponding region is indicated in blue (or 
black, in the case of black-and-white printing). Note that the 
polaronic metallic region with a large z is expressed in red 
(dark gray), while the polaronic phase with smaller z is shown 
in green (light gray). The red (dark gray) and green (light 
gray) regions are continuously connected and the polaronic 
phase with a very small z is defined by the heavy-electron 
state. The bipolaronic insulating phase is expressed by the so- 
lution with an extremely small z and is shown in white in the 
figure. 

The green (light gray) region with a smaller z is shown on 
the left-hand side of the graph, while the white region with 
z — > is shown on the right-hand side of the figure. The blue 



Fig. 8. (Color online) Renormalization factor z on a-/?2 plane for T = 
0.0025. Depending on the magnitude of z, the phase is shown by color gra- 
dation (contrasting density in the case of black-and-white printing) : From 
the left-hand side to the right-hand side of the graph, red (dark gray) for 
0.1 < z < 1, green (light gray) for z on the order of 0.1, blue (black) for 
0.001 < Z < 0.1, and white for z < 0.001 (bipolaronic phase). 



(black) region with the possible heavy-electron state exists be- 
tween them. Note that the narrow blue (black) region with a 
small z on the order of 1/1000 can be found near the bipo- 
laronic state with a > a c2 for large a fa, where a C 2 is the 
upper critical value of the polaron-bipolaron transition. The 
white region cannot be regarded as the heavy-electron state 
owing to the strong localization in the bipolaronic state. 

When the anharmonicity increases with decreasing fa, the 
small-z region expands in the wider region of a around the 
two-solution coexistence or crossover region. When we de- 
crease fa, the expansion of the small-z region is pronounced 
f° r fa ^ —0.25, where the potential shape is changed from 
the fiat single-well type to the double-well type. The present 
results suggest the close relation between the heavy-electron 
state and the anharmonicity. 

3.5 Electron and phonon spectral functions 

Now, we discuss the spectral functions of electrons and 
phonons. For this purpose, we evaluate the DOS obtained 
from the imaginary part of Green's function. The electron 
DOS p e \(uj) and the phonon DOS p p h(u) are respectively 
given by 

—ImG(u + irj), (14) 



and 



Pph(w) 



-ImZ3(w + 177), 



(15) 



where 7/ is a positive infinitesimal quantity, G is defined in 
eq. (10), and the phonon Green's function D is given by 

e -Ei/T _ e ~E s /T 



D(iv n ) 



U\(h + b])\£) 



(16) 



Here, v n — 2irTn is the boson Matsubara frequency. Note 
that we set r\ — 0.004 in order to draw the continuous spec- 
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Fig. 9. Electron DOSs p c i (uj) for three typical potentials as (a) harmonic 
type with /?2 = 0, (b) flat single-well type with 02 = —0.25, and (c) double- 
well type with ft, = -0.5. We set 4 = 0.02 at T = 0.0025. We show the 
a dependence from to 5 from bottom to top for each panel. For each 0.2 
step in the increase in a, the graph is shifted upward by 0.02. 



trum function for clear visibility. 

In Fig. 9, we show the curves for the electron DOS at 
T = 0.0025 for three typical potentials (a) (3 2 = 0, (b) -0.25, 
and (c) —0.5, corresponding to the harmonic, the flat single- 
well, and the double-well types, as shown in Figs. — 1(b)- 1(d), 
respectively. For each panel, we show the curves in the order 
of a from bottom to top. Note that, in the two-solution coex- 
istence region, we plot the results for the polaronic metallic 
solution with a large z. 

At a = 0, a set of two large peaks for high energy uj > 
(low energy uj < 0) and a small peak at uj = can be 
commonly seen in the three panels. Note that the bilaterally 
symmetric graph is due to the particle-hole symmetry. For a 
small a, the structure is almost unchanged. On the other hand, 
for a large a, the set of two large peaks in ui > (a; < 0) 
moves to the right-hand side (left-hand side) of the graph with 
increasing a. Then, the disappearance of the peak at uj = is 
also observed, indicating the disappearance of quasi-particle 
band. The critical values of the disappearance of the peaks at 
oj = are a ~ 4, 2, and 1 for fi 2 = 0, -0.25, and -0.5, 
respectively, which agree well with the a which the transition 
or crossover is indicated. 

Note here that the renormalization factor z discussed above 
should correspond to the decrease in the bandwidth, which is 
well known as the polaronic band narrowing effect. However, 
it is difficult to observe such a band narrowing effect for a 
small a in these graphs, in spite of the significant change in 
z. This is seemingly incomprehensible at a glance, but simi- 
lar results have already been found in the harmonic Holstein 
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OJ 

Fig. 10. Phonon DOSs p p \ l (w) for three typical potentials as (a) harmonic 
type with /?2 = 0, (b) flat single-well type with /?2 = —0.25, and (c) double- 
well type with f3 2 = -0.5. We set /3 4 = 0.02 at T = 0.0025. We show the 
a dependence from to 5 from bottom to top for each panel. For each 0.2 
step increase in a, the graph is shifted upward by 0.04. 



model. 5 1-53 According to previous studies, in the polaronic 
phase, with increasing effect of electron-phonon interaction, 
the central peak at uj = becomes narrower and more pro- 
nounced with the corresponding weight z, while the other 
peaks with total weights of 1 — z exhibit almost no changes. 
The reason why the other peaks do not move at all is that the 
Holstein model has few energy distribution corresponding to 
the upper and lower Hubbard bands, which are split upward 
and downward depending on the Coulomb interaction. Thus, 
the conservation of the DOS for a small a is valid, at least 
within the present calculation on the basis of the exact diago- 
nalization. As for the reproduction of the central peak behav- 
ior, it is necessary to improve the precision of the calculation, 
for instance, by increasing the cluster size. 

Next we move to the phonon spectral function. In Fig. 10, 
we show pph(w) at T = 0.0025 for three typical potentials 
with/3 2 = 0, —0.25, and —0.5. At a = 0, for each three panel, 
distinct peaks can be seen at uj <~ 0.43, 0.3, and 0.04 for f3 2 — 
0, —0.25, and —0.5, respectively. These uj values correspond 
to quasi-harmonic phonons, which agree well with the first 
excitation energies observed in Fig. 2. 

With decreasing (3 2 , the shift of the peaks to the lower- 
energy side is indicated. In the polaronic state for small a, 
the shifts of the low-energy peaks are significant with in- 
creasing a, implying the softening of phonons and the diver- 
gence of charge susceptibility. The heavy-electron tendency 
caused by phonon softening due to anharmonicity has been 
discussed in refs. 54 and 55. In the polaron-bipolaron tran- 
sition or crossover region, the disappearance of the peaks at 
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Fig. 11. (Color online) (a) Lattice fluctuation yj (q 2 ) vs a, (b) renormal- 
ization factor z vs a, (c) double occupancy d vs a, (d) local spin sus- 
ceptibility Xs° c vs a > m d (e) local charge susceptibility Xc° c vs a f° r 
02 = 0, -0.25, -0.5 with T = 0.0025, 0.0035, 0.0045, and 0.0055 
with 04 = 0.02. In the 02 = case, the graph with the crossover result at 
T = 0.0105 is also indicated. Several Xc° c dependences for a large a are 
omitted to avoid overlaps of the graph. Note that xj.° c approaches to a finite 
specific value depending on T for a large a. 



Fig. 12. (Color online) Phase diagrams on a-T plane for (a) 02 = and 
(b) 02 = -0.25 with 0a = 0.02. In the region a c l < a < a c 2, both the 
polaronic and bipolaronic solutions coexist. The dotted-dashed curve shows 
a possible first-order transition temperature and the cross indicates its critical 
point. At lower temperatures, we also plot the expected values of a c i and 
a c 2 with a dotted line. 



lower energy can be seen, suggesting the marked change in 
phonon properties. The critical value of a at which the lowest 
energy peak disappears is in good agreement with the tran- 
sition point or crossover region. In the bipolaronic state, the 
second lower peaks are continuously connected to those in the 
polaronic state and they move to the high-energy side with in- 
creasing a as well as in the electron spectrum cases. 

3.6 Temperature dependence of physical quantities 

Thus far, we have discussed several physical quantities 
focusing our attention on the fixed temperature with T — 
0.0025. In this subsection, we mention the temperature de- 
pendence briefly. Here, we consider three typical potentials 
with /3 2 = 0, -0.25, and -0.5 for /3 4 = 0.02. In Figs. 11, we 
show the a dependences of several physical quantities such 
as (a) lattice fluctuation yj (q 2 ), (b) the renormalization factor 
z, (c) the double occupancy d, (d) the local spin susceptibility 
x' oc , and (e) the local charge susceptibility xi° c for several 
temperatures from T = 0.0015 to 0.0105. 

The marked changes in the physical quantities due to the 
difference in temperature can be seen only in the vicinity of 
the two-solution coexistence or crossover region, and the am- 
plitude of the local charge susceptibility. With increasing tem- 
perature, the reduction in the size of the two-solution coexis- 
tence region is indicated. The edge point for a larger a of the 



two-solution coexistence region a C 2 moves to the left-hand 
side of the graph. As for the shift in a c \, the edge point for 
a smaller a, it is relatively small compared with that in a C 2. 
Furthermore, a c % and a C 2 coincide with each other at high 
temperatures. Then, the two-solution coexistence region dis- 
appears and the first-order transition changes to a crossover. 
Meanwhile, x|° c in the bipolaronic solution for a large a 
rapidly decreases with increasing T. The temperature depen- 
dence of x l ° c m me bipolaronic state exhibits the Curie-law 
behavior x l ° c °c 1 /T, 47 which is similar to that in the case 
with the localized spin in the Mott insulator, where the local 
spin susceptibility shows the Curie law behavior xi° c ozl/T. 
In the harmonic case, the x[° c in the present anharmonic case 
also follows the Curie law. 

In the strong-anharmonicity case with /?2 = —0.5, we point 
out that the crossover behavior is shown at all temperatures, 
even below T = 0.0015 (not shown here). As indicated in 
the lowest-energy phonon spectrum in Fig. 10 and the energy 
levels shown in the bare local phonon potential in Fig. 1(d), 
the low-temperature crossover is attributed to the small exci- 
tation energy of phonons. Namely, the amplitude of the exci- 
tation energy is considered to correspond to the energy scale 
of the polaron-bipolaron transition. Near the crossover region 
at high temperatures, a smooth tail-like behavior with a small 
z is observed. For instance, in the graph for 02 = —0.25, the 
tail-like behavior is observed near a ~ 2.1 for T > 0.0055 
around the crossover region. Moreover, for 02 — —0.5, such 
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Fig. 13. (Color online) First-order polaron-bipolaron transition critical 
point Qcr, T C r vs /?2 for /34 = 0.02. The crosses indicate the temperatures 
where the crossover solution appears on the higher-temperature side and the 
coexistence solutions are found on the lower temperature side. The critical 
points exist on the bars marked by crosses. The expected critical points and 
the projection onto the T = plane are shown as solid and dashed curves, 
respectively. 



behavior can be found at all the temperature, indicating the 
possibility of the heavy-electron state. 

3.7 First-order polaron-bipolaron transition 

In the previous subsection, the two-solution coexistence re- 
gions are indicated for fa = and —0.25 at low temperatures. 
Namely, the first-order polaron-bipolaron transition occurs in 
the region of a c \ < a < a C 2- To discuss this transition, 
we show the a-T phase diagrams for fa = and —0.25 in 
Figs. 12(a) and 12(b), respectively, in which we plot a c i and 
a C 2 for several T obtained from the discontinuity in the lattice 
fluctuation. 

With increasing T, both a c \ and a C 2 show a decreasing 
tendency, with the change in the former being moderate, while 
that in the latter being pronounced. We find that a cl and a C 2 
coincide with each other at a critical point, denoted by a cr 
or T cr . Above this critical temperature, a smooth crossover 
is observed instead of a discontinuous change between the 
polaronic and bipolaronic solutions. 

For fa = and —0.25, the estimated values of T cr are 
~ 0.01 and ~ 0.005, respectively. Below T CT , the first-order 
polaron-bipolaron transition is expected to take place simi- 
larly to that in the case with the Mott transition. 49 Then, we 
plot a possible first-order transition point in the graph. Here, 
we focus on the largeness of the two-solution coexistence re- 
gion, i.e., the region width of \a C 2 — a> c i I an d the critical tem- 
perature T cr height on the a-T plane. It is found that the area 
for fa = —0.25 is smaller than that for fa = in both the a 
width and the T cr height. Since the reduction in the area can be 
associated with the energy scale of the first-order transition, 
this result suggests that the enhancement of the anharmonic- 
ity suppresses the first-order polaron-bipolaron transition. 

Then, we discuss the development of the critical point of 
the first-order polaron-bipolaron transition. In Fig. 13, we 
show the fa dependences of the critical point a cr , T cr . Note 
here that a cr in the crossover region is determined from the 
value at which the curvature of the lattice fluctuation \J (q 2 ) 
changes. In the graph, when fa decreases, both a CY and T cr 



decrease. The decrease in a cr indicates the enhancement of 
the effective electron-phonon interaction coupled to the an- 
harmonicity. As for T cr , we observe a gradual decrease for a 
large fa, while the shape changes to a slightly steeper curve 
for fa ~ —0.25 and finally approaches zero. This T cr be- 
havior seems to be similar to the fa dependence of the first 
excitation energy of phonons in Fig. 2. Thus, the suppression 
effect is thought to be caused by the reduction in the exci- 
tation energy of phonons. For fa < —0.5, we find no evi- 
dence of the two-solution coexistence region for the temper- 
ature region where we have performed our calculations. For 
this reason, we set T cr to in Fig. 13; however, we do not 
exclude the possibility of a polaron-bipolaron transition point 
with an extremely low transition temperature in the region of 
< 

cated in the local phonon problem. 



fa < —0.5, since a small but finite excitation energy is indi- 



4. Discussion and Summary 

In this paper, to discuss the effect of the anharmonicity of 
phonons on the emergence of the heavy-electron state and 
the polaron-bipolaron transition, we have analyzed the an- 
harmonic Holstein model, which describes the interaction be- 
tween conduction electrons and the atom oscillation in the po- 
tential with second- and fourth-order anharmonic terms, using 
the dynamical mean-field theory by the exact diagonalization 
method. We have obtained various physical quantities as func- 
tions of the electron-phonon interaction a for several types of 
anharmonic potential. 

First, we have discussed the effect of anharmonicity on 
the first-order polaron-bipolaron transition. For a weak- 
anharmonicity case with a large fa, with increasing a, we 
have observed an increase in lattice fluctuation, double occu- 
pancy, and local charge susceptibility, but a decrease in the 
renormalization factor z and local spin susceptibility. More- 
over, at low temperatures, evidence of the phase transition 
from the polaronic state to the bipolaronic state, accompanied 
by the changes in the physical quantities, is found for a strong- 
coupling region as the two-solution coexistence region. When 
temperature increases, the region becomes narrower and it fi- 
nally disappears at a certain critical temperature, resulting in 
a smooth crossover between the two states. 

The polaron-bipolaron transition behavior has been re- 
ported in a previous study with harmonic phonons 47 and there 
is no qualitative difference between their tendencies. In a pre- 
vious study with a phonon frequency loq/W = 0.1, the tran- 
sition temperature is estimated to be on the order of 1/1000 
of the bandwidth W. Since the transition is dominated by 
charge fluctuations and it is a relatively spin-independent phe- 
nomenon, we suggest its similarity to the first-order phase 
transition T p = 7.5K in the /3-pyrochlore oxide KOS2O6. 
Here, we note that in this case, the polaronic (bipolaronic) 
state corresponds to the T < T p (T > T p ) side. On the other 
hand, when anharmonicity increases with decreasing fa, we 
have indicated several changes such as the increase in the am- 
plitude of the lattice fluctuation, the decrease in the critical 
value of a, and the transition temperature suppression. Thus, 
in the anharmonic phonon system, a strong electron-phonon 
interaction is not required for the occurrence of the polaron- 
bipolaron transition, even though a lower-temperature envi- 
ronment is necessary. This is considered to be advantageous 
for the explanation of the T p transition. 
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Then, let us discuss how the anharmonicity contributes to 
the heavy-electron state. In this study, we have discussed the 
effective mass of electrons through the renormalization fac- 
tor z as functions of a and fc. With increasing a, z mono- 
tonically decreases, namely, the effective mass increases. The 
heavy-electron state with a small- z solution has been found 
in the narrow region of the polaronic phase near the polaron- 
bipolaron transition. Although smaller values of z have been 
indicated in the bipolaronic solutions, quasi-particles have 
been destroyed owing to the strong localization tendency. 
We have found that the possibility of the heavy-electron for- 
mation is large in the region in which the small-z solu- 
tions are widely distributed. In this sense, the heavy-electron 
state can be found in the region where the crossover is indi- 
cated. The crossover behavior can be seen even in the weak- 
anharmonicity cases with a large /?2 at high temperatures, al- 
though at low temperatures, it becomes a polaron-bipolaron 
transition with a narrow small-z region in the polaronic phase. 
Thus, the crossover region with a high possibility of ex- 
hibiting the heavy-electron state is only found in the strong- 
anharmonicity case with small /3 2 . 

Throughout this paper, we have restricted ourselves to the 
case with the normal state and have not considered the possi- 
ble phase transition to the charge-density-wave (CDW) state, 
which is realized in the Holstein model. Within the DMFT, 
the transition temperature Tcdw to the CDW state in the Hol- 
stein model was found to be 0(1/100) of the bandwidth. 56 As 
the critical temperature T cr of the polaron-bipolaron transition 
obtained in the present study is 0(1/1000) of the bandwidth 
and is much smaller than Tcdw, it seems to be impossible to 
observe the polaron-bipolaron transition in the normal state 
above Tcdw- The effect of the frustration, however, is ex- 
pected to largely suppress Tcdw, while keeping T cr almost 
constant as the polaron-bipolaron transition is exclusively de- 
termined by the local DOS, resulting in Tcdw < T cr . The 
effect of the anharmonicity (3 4 was also found to suppress 
Tcdw except in a small /? 4 regime. 57 In addition, in the 
presence of the Coulomb interaction between electrons, the 
polaron-bipolaron transition in the Hubbard-Holstein model 
was found to be realized for a large-g regime, where T cr 
markedly increases owing to the Coulomb interaction and 
becomes 0(1/100) of the bandwidth, 58 while Tcdw is n °t 
markedly increases but rather reduced. 59 Therefore, we may 
expect that the polaron-bipolaron transition will be observed 
in the normal state above Tcdw in the case with the frustra- 
tion and/or Coulomb interaction. Explicit calculation includ- 
ing the effects of the frustration and the Coulomb interaction 
together with the effect of the doping 60 61 will be an important 
future problem to address. 

In summary, we have investigated the half-filled anhar- 
monic Holstein model using the dynamical mean-field the- 
ory in combination with the exact diagonalization method. 
We have found that, for the weak-anharmonicity case, the 
first-order polaronic -bipolaronic phase transition takes place 
at a critical value of the electron-phonon coupling a, at which 
each physical quantity shows discontinuity. When the anhar- 
monicity is enhanced, the polaron-bipolaron transition tem- 
perature is reduced and the critical value of a decreases. For 
a strong-anharmonicity case, the polaron-bipolaron transition 
eventually changes to a crossover, in which a heavy-electron 
state with a large effective mass is realized owing to the effect 



of anharmonic phonons. 
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